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We present a real-space adaptive-coordinate method, which combines the advantages of the 
finite-difference approach with the accuracy and flexibility of the adaptive coordinate method. The 
discretized Kohn-Sham equations are written in generalized curvilinear coordinates and solved self- 
consistently by means of an iterative approach. The Poisson equation is solved in real space using the 
Multigrid algorithm. We implemented the method on a massively parallel computer, and applied it 
to the calculation of the equilibrium geometry and harmonic vibrational frequencies of the CO2, CO, 
N2 and F2 molecules, yielding excellent agreement with the results of accurate quantum chemistry 
and Local Density Functional calculations. 
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Recently, important efforts in the development of pseu- 
dopotential ab initio electronic structure methods have 
focused on the ability to treat all elements of the peri- 
odic table on an equal footing, and to perform efficiently 
molecular dynamics (MD) simulations. In particular, the 
adaptive coordinate method 1 , proposed to improve the 
efficiency of the widely used plane wave approach, has 
been shown to allow for accurate MD simulations includ- 
ing first-row elements 2 . The applicability of this method 
to the calculation of structural properties of solids 3 and 
to the computation of band structures 4 has also been 
demonstrated. 

Another recent development in electronic structure cal- 
culations, the finite-difference (FD) real-space method 5 , 
has tackled a serious drawback of the plane wave ap- 
proach, i.e. its inability to describe non-periodic systems 
such as clusters or molecules. This problem has been usu- 
ally circumvented by using very large unit cells in order 
to minimize the interaction between replicas of the sys- 
tem. However, this often leads to very slow convergence 
of the results as the cell size is increased, and can be for- 
bidding for charged systems. The FD real-space method 
is based on a discretization of the Kohn-Sham equations, 
and does not imply any assumption of periodicity of the 
solutions, thus allowing one to choose boundary condi- 
tions (periodic or Dirichlet) most appropriate to the sys- 
tem under study. It was pointed out in Ref. 5 that the 
finite-difference approach might be well adapted to mod- 
ern parallel computer architectures since it avoids com- 
pletely the use of Fourier transforms. We note that an 
important feature of real-space methods is that they can 
lead to a straightforward implementation of linear scaling 
electronic structure algorithms (0(N) methods) within a 
self-consistent Density Functional approach 6 . Such algo- 
rithms can be naturally formulated in real space, since 
they make use of localized electronic orbitals defined in 
domains smaller than the entire simulation cell. 

For the treatment of first row elements, the finite- 
difference approach suffers from the same inefficiency 
problem as plane waves. Indeed, in order to system- 
atically improve the accuracy of a calculation, the grid 
spacing h must be reduced on the entire grid, even though 



an increased resolution would be needed only in regions 
of rapidly varying potential. In the case of first-row 
elements this can be computationally very demanding. 
However, a local reduction of the grid spacing can be 
achieved by applying the concept of adaptive curvilinear 
coordinates to the finite-difference approach. 

In this paper, we present a real-space adaptive- 
coordinate method, which combines the advantages of 
the FD approach with the accuracy and flexibility of 
adaptive coordinates. The discretized Kohn-Sham equa- 
tions are rewritten in generalized curvilinear coordinates 
and solved self-consistently using an iterative approach. 
The method is implemented on a massively parallel com- 
puter, and applied to the calculation of the equilibrium 
geometry and harmonic vibrational frequencies of the 
CO2, CO, N2 and F2 molecules, yielding excellent agree- 
ment with the results of other accurate first-principles 
calculations. Our results show that the use of curvilin- 
ear coordinates is essential to obtain a sufficiently small 
grid spacing in the vicinity of atoms, while keeping the 
numerical effort limited. 

Curvilinear coordinates can be chosen to be either fully 
adaptive 1,3,4 , or defined from given deformation func- 
tions as in Ref. 2, which is more efficient in molecular 
dynamics simulations. We follow Ref. 2, and define a 
coordinate transformation £ — ► x(£) by 



? =X l +J2(x l -K)f a (\x-Ra 



(1) 



where R a denote the atomic positions and the deforma- 
tion functions f a (r) are given by 



fa(r) = A a — tanh — exp 

r a„ 



(2) 



This coordinate transformation maps a regular grid in 
^-space onto a curvilinear grid in (real) x-space, and re- 
duces the grid spacing in x-space near atoms. In Eq. 2 
a a is the range of enhancement of the resolution around 
atom a, and b a determines a distance over which Eu- 
clidean coordinates are recovered. The local reduction in 
grid spacing at the center of an atom of species a is given 
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by h e ff/ho = 1/(1 + A a ), where ho denotes the grid spac- 
ing in Euclidean coordinates. The corresponding local 
effective energy cutoff at the atomic site is 
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The wavefunctions are then written as the product 7 

dc 1 



t/>(x) = <f>(-x) 



dx 



(4) 



where \d£/dx\ is the Jacobian of the transformation 
x — ► £. The function <f>, which was expanded in Fourier 
series in the plane-wave adaptive coordinate method 1 is 
now represented by its values on the curvilinear grid. The 
product representation Eq. (4) of ip has the advantage of 
conserving scalar products during a change of coordinate 
transformation x(£) —> x'(£). This has the important 
consequence that a set of orthonormal wavefunctions re- 
mains orthonormal when x(£) changes due to the motion 
of atoms in a molecular dynamics simulation. 
We write the curvilinear Laplacian as 



a ® d 



d 2 C d 



d£, 1 d£, 3 dx k dx k d£, 1 
where we expand the second term on the r.h.s. as 
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dx k dx k dx k d^ 1 \dx k 



(5) 



(6) 



Summation over repeated indices is implied throughout 
this paper. The discretized version of A is obtained by 
replacing all derivatives with respect to !;* in (5) and (6) 
by finite-difference formulae and by using the analytic 
expression of d^' /dx 3 derived from the definition of £(x). 
In Eq.(5), g' 3 is the inverse of the Riemann metric ten- 
sor and is obtained directly from the analytic expression 
of d^ l /dx 3 . Note that in Euclidean coordinates, the dis- 
cretized Laplacian is self-adjoint and positive-definite to 
all orders in h. On the contrary, the discretized curvi- 
linear Laplacian is only self-adjoint and positive-definite 
up to terms of the order of the truncation error (i.e. 
0(h N ), where N is the order of the finite-difference for- 
mula). In practical applications, this departure from self- 
adjointness and positive-definiteness has no noticeable 
effect on the convergence of the results, if the deforma- 
tion functions /„ are sufficiently smooth. The trunca- 
tion error can however become significant if the Jacobian 
d^' /dx 3 is varying too rapidly on the scale of the grid 
spacing. We use finite-difference formulae of order 2 and 
4 for the discretized derivatives i.e. 3- and 5-point for- 
mulae in each direction. Higher order finite difference 
formulae are not needed since the reduction in grid spac- 
ing provided by the adaptive coordinates usually brings 
sufficient accuracy. Low-order formulae also reduce the 
amount of interprocessor communications on a parallel 
computer. 



The potential energy operator consists of norm- 
conserving, non-local pseudopotentials in the Kleinman- 
Bylander 8 form. The integrals associated with the local 
potential, non-local potential, and exchange correlation 
energy are evaluated by straightforward sums over the 
grid points, appropriately weighted with the Jacobian. 

The properties of the curvlinear discretized Hamilto- 
nian differ in many respects from those of the plane- 
wave Hamiltonian. First, the wavefunctions, the charge 
density and the potentials are all represented on the 
same grid, whereas in the plane-wave method, the po- 
tentials and the wave-functions are usually represented 
as Fourier series having different energy cutoffs. Fur- 
thermore, finite-difference methods do not rely on an ex- 
pansion of the solutions on a basis set, and the energy 
therefore does not obey a variational principle. Instead, 
the truncation error of the energy can have an arbitrary 
sign, and in practice, the convergence of the energy to the 
exact value as h —> is often attained from below. Fi- 
nally, unlike in the plane-wave approach, the discretized 
Hamiltonian is not in general invariant under uniform 
translations of the system with respect to the position 
of the grid. This loss of translational invariance requires 
special care in the convergence of the calculation if spu- 
rious forces are to be avoided. In Euclidean coordinates, 
it is easily shown 9 that translational invariance of the 
energy is restored if the ionic pseudopotential is band- 
width limited, i.e. if its Fourier transform V(q) vanishes 
beyond the Nyquist critical wave vector q c = ir/h. In 
practice, since typical norm-conserving pseudopotentials 
are not bandwidth limited, translational invariance of the 
energy can be recovered by either reducing h until V(q c ) 
becomes negligible, or by appropriately filtering the po- 
tential. Although no simple equivalent to the Nyquist 
critical wave vector q c exists in curvilinear coordinates, 
we found that a local reduction of the grid spacing h 
near atoms, i.e. where the ionic potential has rapid os- 
cillations, can be sufficient to restore the translational 
invariance of the energy. 

Within the adaptive coordinate approach, the Poisson 
equation must be solved by an iterative method, since 
the curvilinear Laplacian is not diagonal in either real 
or reciprocal space. We solve Poisson's equation entirely 
in real space where the discretized curvilinear Laplacian 
is sparse. This allows one to work in domains of ar- 
bitrary shapes, with arbitrary boundary conditions. It 
also avoids the use of Fourier transforms, which is ad- 
vantageous on massively parallel computers. As an itera- 
tive procedure, we choose the Multigrid method 10 , which 
has the property of including an implicit precondition- 
ing on all length scales, and is therefore rapidly conver- 
gent. In order to ensure the convergence of the Multi- 
grid iteration, the discretized Laplacian operator must 
be non-singular. This can be accomplished by project- 
ing appropriately the solution out of the nullspace of the 
discretized Laplacian. For our choice of the discretized 
Laplacian Eq.(5), and for periodic boundary conditions, 
the nullspace reduces to the constant function. 



2 




2.05 2.1 2.15 2.2 2.25 2.3 2.35 

C-0 distance (a.u.) 

FIG. 1. Total energy of the CO2 molecule as a function of 
the C-O distance in the linear conformation, calculated using 
Euclidean coordinates (open circles) and adaptive curvilin- 
ear coordinates (dots). Spline interpolations are shown as a 
guide to the eye. The two curves are shifted arbitrarily in 
order to allow for their comparison. The experimental value 
of the equilibrium C-O distance (arrow) and the grid spacing 
in Euclidean coordinates are also indicated. 

The implementation of the real-space adaptive- 
coordinate approach was carried out on a Cray T3D mas- 
sively parallel computer, using the Parallel Virtual Ma- 
chine (PVM) environment. For simplicity, the simulation 
cell was divided into domains of equal shape and size 
which were distributed among processors. Other choices 
could be made to improve e.g. load balancing of the com- 
putations. Communications are limited to nearest neigh- 
bouring processors, and only imply the transfer of data 
defined on domain boundary layers. In particular, no 
global communication of a function defined on the entire 
grid is ever needed. 

We applied the real-space adaptive-coordinate ap- 
proach to the calculation of the equilibrium configura- 
tion and of the harmonic vibrational frequencies of the 
CO2, CO, N2 and F2 molecules, within the local density 
approximation (LDA) of density functional theory. The 
experimental 11,12 equilibrium bond lengths in CO2, CO 
and N2 are 2.19, 2.13 and 2.07 a.u. respectively. Repro- 
ducing such short bond lengths, as well as vibrational 
frequencies in these molecules represents a stringent test 
for pseudopotentials. The description of F2 is also a diffi- 
cult test for a pseudopotential approach since the fluorine 
p wavefunction is one of the most strongly localized va- 
lence wavefunction of all elements in the Periodic Table. 

For these reasons, we used the norm-conserving pseu- 
dopotentials of Bachelet, Hamann and Schliiter 13 , which 
reproduce accurately the atomic all-electron valence 
wavefunctions already at distances of 0.8, 0.7, 0.6 and 
0.5 a.u. from the nuclei, for carbon, nitrogen, oxygen 



and fluorine respectively. We first computed the total 
energy of CO2 as a function of the C-O distance c?co 
using Euclidean coordinates and a grid of 96 x 72 x 72 
points in a cell of 20 x 15 x 15 a.u. This resolution would 
correspond to an effective energy cutoff of 227 Ry in the 
expansion of wavefunctions, charge density and poten- 
tials in a conventional plane wave calculation. We used a 
5-point finite-difference formula in each direction to com- 
pute derivatives. The resulting total energies are shown 
on Fig. 1 (open circles), together with a spline interpola- 
tion passing through calculated points (solid line). The 
loss of translational invariance manifests itself strikingly 
by the presence of two minima in E(dco) in the vicinity 
of the experimental value of c?co- The distance between 
the two minima is close to the grid spacing ho, which con- 
firms that the oscillations in E(dco) are spurious and de- 
pend on the relative positions of the atoms with respect 
to grid points. This clearly demonstrates that the res- 
olution used in Euclidean coordinates is insufficient for 
a calculation of the equilibrium value of c?co- We then 
repeated the calculation using the same number of grid 
points and the same unit cell, but with curvilinear coor- 
dinates, defined by the deformations given in Eq. (2) with 
Ac = 0.258, ac = 1.10 a.u. and be = 4 a.u. for carbon, 
and Aq = 0.989, ao = 0.60 a.u. and 6 = 4 a.u. for oxy- 
gen. This corresponds to effective energy cutoffs on the 
atomic sites of 360 Ry and 900 Ry for carbon and oxygen 
respectively. The total energies obtained with curvilin- 
ear coordinates are also shown in Fig. 1 (dots), together 
with a spline interpolation (solid line). The curve E(dco) 
has only one minimum at c?co = 2.198 a.u. which com- 
pares nicely with the experimental value of the equilib- 
rium bond length of 2.192 a.u. The harmonic vibrational 
frequency of the <r+ symmetric stretching mode was ex- 
tracted from a fourth order polynomial fit to E(dco)- A 
similar calculation was carried out to determine the har- 
monic frequencies of the ir u bending mode and of the <r+ 
asymmetric stretching mode. We note that harmonic fre- 
quencies cannot be compared directly with experimental 
frequencies, due to the presence of anharmonic terms in 
the potential energy surface, especially for the <r+ mode 
which is split by a Fermi resonance 14 . They can be com- 
pared, however, to harmonic frequencies obtained from 
a parametrization of the potential energy surface of CO2 
adjusted to experimental spectra 11,15 . Our calculated 
frequencies are 1336, 648 and 2374 cm -1 for the <r+, ir u 
and <7+ modes, respectively, which compare well with the 
corresponding values extracted from experiment (1354, 
673 and 2397 cm -1 ). A summary of the results for CO2, 
and a comparison with other calculations are given in 
Table I. 

We also computed the equilibrium bond lengths and 
the harmonic vibrational frequencies of the CO, N2, and 
F2 molecules. For N2, and F2, we used deformation func- 
tions /„ yielding an effective energy cutoff of 800 and 990 
Ry at the atomic sites respectively. 
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TABLE I. Equilibrium C-O distance dco and har- 
monic vibrational frequencies in CO2 calculated using the 
real-space adaptive-coordinate method, compared with har- 
monic frequencies adjusted to experimental spectra, and 
with the results of other calculations (LDA, Hartree-Fock 
(SCF), configuration interaction (CISD) and coupled cluster 
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This work 
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1336 
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2.192 
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a From Ref. 17. 
b From Ref. 11. 



Our results are summarized in Table II where they are 
compared with those of other LDA calculations 16,17 and 
with experiment. The overall agreement with other LDA 
results is excellent. The agreement with experiment is 
also remarkably good, except for F2, for which a more 
accurate treatment of correlation is needed. 

In conclusion, we have presented a real-space adaptive- 
coordinate method for electronic structure calculations. 
This approach combines the advantages of the finite- 
difference method, with the capability of the adaptive- 
coordinate method to treat first-row and heavier elements 
on the same footing. We have shown that the Poisson 
equation in curvilinear coordinates can be solved effi- 
ciently using the Multigrid method, and that a local re- 
duction of the grid spacing near atoms can restore the 
translational invariance of the energy. The method was 
implemented on a massively parallel computer, and ap- 
plied to the calculation of the equilibrium bond lengths 
and vibrational frequencies of several molecules contain- 
ing first-row elements, yielding excellent agreement with 

TABLE II. Equilibrium bond lengths d and harmonic vi- 
brational frequencies co e in CO, N2 and F2 calculated using the 
real-space adaptive-coordinate method, compared with other 
LDA results and with values extracted from experiment. 
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2387 a 
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This work 


2151 


2382 


1051 


Expt. 


2170 c 


2358 a 


892 a 



a From Ref. 16. 
b From Ref. 17. 
c From Ref. 12. 



other accurate first principles results. Work is in progress 
to apply the method to ab-initio molecular dynamics sim- 
ulations, and to extend it to a linear scaling formulation. 
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